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Abstract. In this paper we examine the reconstruction of functions of any dimension 
from hyperplanar projections. This is a generalization of a problem that has 
generated much interest recently, especially in the field of medical imaging. 
Computed Axial Tomography (CAT) and Nuclear Magnetic Resonance (NMR) are 
two medical techniques that fall in this framework. CAT scans measure the x-ray 
density along lines through the body, while NMR scans measure the hydrogen 
density along planes through the body. 

Here we will examine reconstruction methods that involve backprojecting the 
projection data and summing this over the entire region of interest. There are two 
methods for doing this. One method is to filter the projection data first, and then 
backproject this filtered data and sum over all projection directions. The other 
method is to backproject and sum the projection data first, and then filter. The 
two methods are mathematically equivalent, producing very similar equations. 

We will derive the reconstruction formulas for both methods for any number of 
dimensions. We will examine the cases of two and three dimensions, since these are 
the only ones encountered in practice. The equations are very different for these 
cases. In general, the equations are very different for even and odd dimensionality. 
We will discuss why this is so, and show that the equations for even and odd 
dimensionality are related by the Hilbert Transform. 
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1. Introduction 


In this paper we will examine the reconstruction of functions of any dimension 
from hyperplanar projections. This is a generalization of a problem that has 
generated much interest recently, especially in the field of medical imaging. 
Computed Axial Tomography (CAT) and Nuclear Magnetic Resonance (NMR) are 
two medical techniques that fall in this framework. 

CAT scans measure the x-ray density along lines through the body. Reconstruc¬ 
tion of the density distribution results in a 2-dimensional image from line projections. 
The line projections are degenerate hyperplanar projections for the 2-dimensional 
case. 

NMR scans measure the hydrogen density along planes through the body. 
Reconstruction results in a 3-dimensional density distribution. The reconstructed 
distribution is usually sliced along an arbitrary plane to produce a 2-dimensional 
image for display purposes. 

Here we will examine reconstruction methods that involve backprojecting the 
projection data and summing this over the entire region of interest. Such methods 
produce a smeared image, so it is necessary to undo this smearing. There are 
two methods for doing this. One method, known as the rho-filtered layergram, 
is to filter the projection data first, and then backproject this filtered data and 
sum over all projection directions. The other approach is to backproject and sum 
the projection data first, and then filter. The two methods are mathematically 
equivalent, producing very similar equations. They differ in that in the first case 
the filtering need only be performed in one dimension, while in the second case, 
filtering must be performed in N dimensions. One dimensional filtering is generally 
easier to implement. 

We will derive the reconstruction formulas for both methods for any number 
of dimensions. We will examine the cases of two and three dimensions, since these 
are the only ones encountered in practice. The equations are very different for these 
cases. In general, the equations are very different for even and odd dimensionality. 
We will discuss why this is so, and show that the equations for even and odd 
dimensionality are related by the Hilbert Transform. 
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/“\ 2. Convolution Followed by Backprojectionand Summation 

Following (Louis and Natterer 1983), we first introduce some notation. Let R N 
be the N-dimensional real space, and let S N ~ l be the set of directions in R n . s n - 1 
is formed from all unit vectors in R N . 

S N ~* = {x £ R N such that jx| = 1} (2.1) 


Let us consider the following problem. The density of the object under study 
will be denoted by /(x), x £ R N . We are given projections p(s,n), n £ S N ~' and 
wish to reconstruct /(x). Figure 1 illustrates the geometry of the problem for CAT 
and NMR. p(s, n) is defined by 


p(s,n)= / /(x) dx 

s=xn 

= f /(x)6(s — x • n) dx (2.2) 

R” 




The Fourier Transform of /(x) and its inverse are given by 

F{u) = f /(x)e- J ' wx dx (2.3) 

R N 

/(x) = (2 tt)~ n f F{oj)e 3U,x ckj (2.4) 

R N 


The problem of reconstruction from hyperplanar projections is most naturally 
handled in polar coordinates because the arguments to the projection p(s, n) can 
be thought of as polar coordinates in projection space. In N-dimensional polar 
coordinates, the Fourier Transform of /(r, m) and its inverse are given by 

F(p,a) — 2 -1 j I f(r >m )e-^ a m \r\ N - l drdm (2.5) 

S N ~ 1 R 

f(r, m) = (2tt)- 7V 2- 1 J j F{p,a)e jr f ,a - m \p\ N - 1 dp da (2.6) 

S *- 1 R 

The N-dimensional Fourier Transform of /(r, m) and the one-dimensional 
Fourier Transform of the projection p(s,n) are related by the projection theorem 
(Appendix A). 

F(p, a) = P(p, a) (2.7) 


2 
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rs 


So, if we examine the inverse Fourier Transform of /(x) in polar coordinates, 
we get 


/(x) = ( 2tt)- n 2- 1 I I F(p, a )e^ ax \p\ N - 1 dpd, 
S N ~ l R 

J \(2ir)- N 2~ 1 J P{p,c)e>'’ a ^\p\ N - 1 dp 


S N-l 


da 


(2.8) 


We see from (2.8) that /(x) can be obtained by computing the expression in 
brackets, and smearing this function over all directions. This smearing is simply the 
backprojection-summation operation. Most authors (e.g. Chiu et.al. 1980, Lewitt 
1983, Shepp 1980) normalize the backprojection operation by dividing by the 
surface area of a unit sphere in N dimensions. Let An be this factor. 


An= J 

S N-1 


„ N 
, 2?rT 

da — - 

r(f) 


(2.9) 


With the inclusion of this normalization term, we can rewrite /(x) as 







1 

An 


/ 

s N - 1 


A n (2'k)~ n 2~ 1 


J P(p,a)e ipax 

R 


\P\ N ld P 


da 


( 2 . 10 ) 


The bracketed expression is the inverse FT of a product. It can also be expressed 
as a convolution. 

A N (2ir)~ N 2~' / P(p, a ) c >^\p\ N -'-dp 
R 

= (27T)- 1 / [T(f )]“ 1 27rTP(p, a )(27r) 1 - A '2- 1 |p|^- 1 e J > a - x dp 
R 

— p(a • x,a) 0 FT'“ 1 |[r(^)]~ 1 2 1_Ar 7 r 1_ T|p| Ar “ 1 | where 0 = convolution 
= p(s,a)0 ff (s) (2.11) 

Where g(s) = Fr- 1 |[r(f)]~ 1 2 1 - iV 7r 1 -^|p| ;v - 1 J (2.12) 

with Fourier Transform G(p) — [T(y )] _1 2 1 ~ A( 7r 1 ~'2"|p| J?v '~ 1 (2.13) 

We will treat the cases of N even and odd separately. 

Case N Even 

The fact that N is even can be used to advantage by introducing the cosine 
transform. 


4 
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/ °o 

-00 

poo 

= [r(frV^TT-f / p N ~\e jsp + e~ jsp )dp 

JQ 

p oo 

= [r(fr'a'-'V-f / o p "- 1 cos(v)^ 

= [r(f)] _1 2 1 - Ar 7r-Tr(Ar)cos(^)s _iV 

= 2 1 - JV 7T-^^^(-l)^s- JV (2.14) 

1 l 2 i 


Strictly speaking, the integral only converges for 0 < AT < 1 (Gradshteyn and 
Ryzhik 1980). It has been argued using convergence functions (Horn 1973) that the 
same form holds for all N. This is not quite correct. 

This can be seen by considering the total area of g(s), f^x, g(s) ds. This must 
equal 0 since G(p = 0) = 0 for N even and positive. But, when N is even, 1 /a N is 
always positive and g(s) cannot possibly integrate to 0. 

This discrepancy can be resolved by taking g(s) to be a generalized function 
(Gel’fand and Shilov 1964, Lighthill 1958) rather than an ordinary function. Thus, 
g(s) can be defined as 


Where g e {s) 


g i s ) = Jim g c {s) 


if |s| > 6 
if Is I < e 


(2.15) 


(1-A0« N 

To simplify notation, we will write s~ N to denote the generalized function 


lim ( ! 

e ~*° [ (1 -N)e" 


if |s| > e 
if |s| < e 


(2.16) 


We will continue to write ^ for the ordinary function. 

The deblurring function blows up rapidly at the origin, for all values of N. 
This may introduce numerical problems when trying to actually compute the 
convolution. 

Notice that p(s) has an infinite region of support. Therefore the reconstruction 
of /(x) at a single point will depend on all values of the projection p(s,n). 
Reconstruction in this case is not local. 

The above formula is only valid for N even. When N is odd, formally applying 
the formula yields 0 everywhere, with the possible exception of the origin, because 
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the factor of cos (^r) is 0. At the origin, we get 0 times a function that goes up 
as s~ N . Therefore the case of N odd must be treated separately to get meaningful 
results. 1 

Case N Odd 

When N is odd we have from (2.13) 

G{p) = [r(f )]- 1 2 1 - jV 7T 1 “T^- 1 (2.17) 


The differentiation operator has transform jp, and the N — 1 derivative, N 
odd, has transform p N ~ 1 (— 1) * • So, 


g{s) = 


r(f) 


iV-1 

~T~ 


d N ~' 
ds N ~ 1 


(2.18) 


In this case t/(s) is a differential operator. Its region of support is limited. 
Therefore the reconstruction of /(x) at a single point will depend only on values of 
the projection p(s, n) that include the point of interest. Reconstruction in this case 
is local. 

Summary 

In the convolution-backprojection method, each projection must be filtered 
before backprojection and summation. The filtering function is 


g(s) = 
In CAT, N = 2 

In NMR, N = 3 


f2 1 -^-fpl(-l)f S ^ 


,1-TV^l-f _1_ ' 


r(f) 


(-i Y 




1 -2 


?(S) = _ 27T S 


?(s) 


27r ds 2 


if N even 
if N odd 


(2.19) 


( 2 . 20 ) 


( 2 . 21 ) 


These special cases agree with (Chiu et.al. 1980), (Shepp 1980), and 
(Tanakaffil979). 


Actually, if we consider g(a) to be a function of N, it happens that it is an entire function, 
and therefore can be defined even when N is odd. Such an approach is taken in (Gel’fand and 
Shilov 1964) to obtain a differential operator when N is odd. 
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3. Backprojectionand Summation Followed by Convolution 

In this method, we first form the backprojection of the projection data, sum 
over all directions, and then filter. As before, the backprojection is normalized by 
A]y. The backprojection-summation b(x) is obtained by 


6(x) = A^ 1 / p(x • n, n) dn 

S S-1 

— j J f(x')6(x • n — x' • n) dx'dn 

S N-1 RN 

J /(x')Aj^ 1 J 6(x ■ n —x' ■ n)dn dx! 


RN gN~ 1 

/(x) ® h(x) 


(3.1) 


The backprojection-summation operation is equivalent to convolving the density 
function /(x) with the blurring function h(x). We need to find what h(x) is. In 
Appendix D it is shown that 


[ <5(x • n) dn 

S N ~ 1 

Therefore the blurring function is 

'N 


_ AT-1 

27r~r- 


|x|r(^V) 


N-l 


(3.2) 


27T 2 
1 


Mx, - 

" ~ WUV) 
_ r(f) 

x| \/i : r(^ri) 


(3.3) 


We would like to find another function, g(x) that will undo the blurring caused 
by h(x). g must be the convolutional inverse of h i.e. g(x)<$Z) h(x) = <5(x). We observe 
that since h is symmetric, g must also be spherically symmetric. Therefore, in this 
section we may take g — g(r). The best way to find g is to take Fourier Transforms. 
From above, we must have 

G ^=m (3 - 4) 

Using the result in Appendix C concerning the FT of |r|, we have, with 
k==-1 




GW 


2' v ->7r*-‘r(£)H 


l-N 


lN-1 


r(f) 


(3.5) 

(3.6) 


As before, we will treat the cases of N even and odd separately. 
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/■s 


/""■N 




Case N Even 

In the case of N even, we can use the results of Appendix C to determine the 
inverse FT of G(p). 


g(r) = 2'- N r l -t^—2 

r (f) 
i_ N r(;v-j) 
~ ' r(f)r(i^ 


N -1 


7T 


tv N-i-N 
—2- 


r( 


-N+U 
2 1 


.1-2JV 


r -N-A^+l 


(3.7) 


Most of the comments that were made about convolution-backprojection for 
N even also apply here. Specifically, the deblurring function blows up at the 
origin. Also, its region of support is infinite rather than local. As before, g(r ) is a 
generalized function. 

The above formula is only valid for N even. When N is odd, the factor r(^-y^) 
is undefined. Therefore, the case of N odd must be treated separately. 

Case N Odd 

This is solved in exactly the same manner as in the convolution-backprojection 
method. When N is odd we have from (3.6) 

G(p) — 

Uf) 


(3.8) 


The Laplacian operator, V 2 , has transform — p 2 . So, 

gtr) = 2 -5— (-1)^ V w “‘ (3.9) 

r(£) 


Where 1 indicates to take the Laplacian times. This is equivalent to 

x2k 

v 2fc = EE---E--- (3.io) 

frt tdxidxi'.dxi 

This is the same form as was obtained with the convolution-backprojection 
method, except the differentiation in one dimension has been replaced by a series of 
Laplacians. As before, g(r ) is a differential operator. Its region of support is limited. 
Therefore the reconstruction of /(x) at a single point will depend only on values 
of the backprojection 6(x) in a small neighborhood about the point, which in turn 
depend depend only on values of the projection p(s,n) that include the point of 
interest. Reconstruction in this case is local. 


8 
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Summary 


In the backprojection-convolution method, each projection must be filtered 
after backprojection and summation. The filtering function is 

g{r) = < 

f -1-JV r (^-|) -1-2JV 

2 1 ~ A V-f 

if N even 

if N odd 

(3.11) 

In CAT, N = 2 

s'" 


(3.12) 

In NMR, N = 3 

9 (X) = -^ 


(3.13) 


These special cases also agree with (Chiu et.al. 1980), (Shepp 1980), and 
(Tanaka 1979). 
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/""s 4. Comparison of Even and Odd Dimensional Reconstruction 

In this chapter we will discuss why the reconstruction equations look so 
different in even and odd dimension. To make the discussion concrete, consider the 
convolution-backprojection method. To repeat, we have (2.13) 

r(f) 






The filtering function is (2.19) 


9{s) 


f2^-ff{fl(-l)f S -^ 




AT 
‘ 2 


fc(-l) 


N — l iN~ 1 


r(f) 


ds N ~ l 


if N even 
if N odd 


It is possible to explain the differences by noting that 1 = p N ~ i when 
N is odd. This leads to the second form above, in which the filtering operation is 
equivalent to differentiation. When N is even, then sgn(p). Since 

multiplication in the transform domain corresponds to convolution in the spatial 
domain, we can try to perform the convolution directly. So let us redo the analysis 
for N even, using this approach. 

For N even, we have 


CM = 2 1 --'V-¥^-/ r - 1 sgnM 

r(f) 


Now the inverse transform of p N 1 is 


ft - V 1 } = (- y )^- 1 


d*- 1 

ds^* 1 


(4.1) 


(4.2) 


And the inverse transform of sgn(p) is 


WHsgnM) 


ITS 


Putting it all together, we get 


ff(«) 


,1-V^l-f. 1 ( -\N -1 dN 1 


r (f) 


-(-/) 


ds"-' 
N-l-N 


2 1 - n *-2 ^(N - 1)! a- 

N' 


ns 
N -1 


■N 


r (f) 


,1 —M _AfT(iV) / .AT _AT 

* 2 !) 2 s 


m) 


(4.3) 


(4.4) 
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This is exactly the result obtained through direct transformation using the 
formula in Appendix C. This shows that the basic difference between even and odd 
dimensions is that raising the magnitude of p to an odd power is equivalent to 
differentiating the function £ in the spatial domain. To make the similarity more 
apparent, we can rewrite the correction function as 


ff( s ) = { 


# 


__L_ »AT—1 
T(f) J 
l yiV-l 

r(f) ; 


d N ~ 1 

ds"- 1 

d"- 1 



if N even 
if N odd 


(4.5) 


The above expressions for g[s ) form a Hilbert Transform (Bracewell 1965) pair. 
The Hilbert Transform of an even function is an odd function, and the transform of 
an odd function is even. Thus the Hilbert Transform can be regarded as correcting 
the antisymmetric differentiation when N is even. 

This analysis does not carry over to the backprojection-convolution case. This 
is because it is not possible to define sgn(p) in more than one dimension. Nor can 
the Hilbert Transform be defined. However, the underlying principle is the same. 
That is, the correction function must be spherically symmetric in all cases. When 
the dimensionality of the problem is even, this leads naturally to derivatives of 
even order. When the dimensionality of the problem is odd, we are not allowed to 
use the odd derivatives, because they are not symmetric functions. 


Acknowledgment: I are grateful to B.K.P. Horn for suggesting this field of 
research, and for his many valuable suggestions and discussions. Thanks also to A. 
Yuille for his comments and suggestions. 
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Appendix A 

In Appendix A we prove the projection theorem (Mersereau and Oppenheim 
1974) in polar coordinates. Let p(s,a) be the projection of /(x) along a plane 
perpendicular to a. Then the one-dimensional Fourier Transform of p and the 
N-dimensional Transform of / are related as follows. 

P(p,a) = J p{s,a)e~ 3Sp ds 

R 

— J J f(x)8(s — x-a)dxe~ JS(, ds 

R R» 

— f 2 _1 / f f(r,m)8(s — rm ■ a)\r\ N ~ 1 dr dme~ 3 ' sp ds 

R S N ~ 1 R N 

— 2 _1 J / /(r,m)|r| Ar ~ 1 e” jrpm ' Q! drdm 

S N -1 R 

— F(p, a) (A.1) 

The FT of p(s,a) equals the FT of /(x) evaluated in the a direction. 
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Appendix B 

In Appendix B we derive the Fourier Transform of a spherically symmetric 
function and the inverse transform. Let /(x) = f(r) with Fourier Transform 
F(u )) = F(p). In polar coordinates the elements of x and w can be defined by 

xq — r sin $o- • • sin 2 
N-2 

Xk — r cos dk-i JJ sin 1 < k < N — 2 
i—k 

XjV-l — PCOS0/V-2 




u >0 — p sin q:o' • • sin ctN-2 
N-2 

Uk = pcosotk-i JJ sin a,-, 1 < k < N — 2 

i=k 

WJV-l = P COS OL jV—2 

The Fourier Transform is given by 

F[ui)= f fMe-^dx 

R N 


(B.l) 


-jrp 

/ OO /*7T /*7T 

-oo/o -■/„ 


/sin 5o - " sin sin ao --- sin aw _2 + \ 

cos 9q sin &i■■■ sin 8^-2 c°s ao sin «i — sin cxn -2 
4- ••• + 

\cos cos ajv-2 / 


r N 1 sin 1 Oi-• ■ sin^ 2 6n-2^0' • dr 




(B.2) 

The transform of a spherically symmetric function is also spherically symmetric. 
So we can take a,- = 0 without loss of generality. The integration over r need only 
be performed from 0 on up, as long as a factor of 2 is included. Using the facts 
that 

sm k 0(16 = ^ 

m + i) 

nt) = v* 

and 
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We have 


r°° r v t* fir 

F (P) = 2 J 0 f{ r > N ~ X dr] Q sm o 0ode Q .-.J Q sin*" 3 9 N _ 3 dO N ^ J e~^ cos6 ^ s in N ~ 2 9 N ^ 2 d9 N - 2 


r oo ;v —6 f 

2 Jo f( r y N ~ ldr n (** 

i=o V 

/* QQ 

(27r)fp 1_ f / r%J&_ 1 (rp)f(r)dr 




(B.3) 


Interchanging r and p and dividing by {2tc) n gives us the inverse Fourier 
Transform 


p oo 

/(r) = (2ir)“T r l ~! J pf J^_ l (rp)F(p) dp 


(B.4) 






14 
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,/«*s 




Appendix C 


In Appendix C we apply the spherically symmetric form of the Fourier 
Transform to the function f[r) = r k . Since r — |x| and p = |o»|, we will assume in 
this appendix that r and p are non-negative. From Appendix B we have 


Let z — rp 


*00 

F (p) ~ (2n) r p 1 ~% J o JN_ x {rp)f{r) dr 

N hi hi 

{2-k)~z p l ~~* / r~z Jri_ 1 (rp)r k dr 

JQ 2 

* oo 

( 2 tt ) t P~ k ~ N f Q z% +k Jz^z) dz 


2 N+k -N-k 

n-i) 


Using the fact that 


(C.l) 


F(p) 


fOO TV P+g+l \ 

/ x q J p [x) dx — 2 ? —-——p 

JO pv ' p( P-?+ 1 ) 

Interchanging r and p and dividing by (2^)^ gives us the inverse transform of 


2 > 


r(-|) 


(0.2) 




15 





Gennert 


Reconstruction from Projections 


Appendix D 


In Appendix D we derive an expression for J s n-i 6{pc • n) dn This integral is 
obviously independent of the direction of x, but is dependent on its magnitude. 
We can express n in polar coordinates, and let x lie along the last component of n, 
since the orientation of n is arbitrary. Define n by 

no = sin 9q- • • sin 0/v-2 
N-2 

nk — cos 0fc_i JJ sin 9{, 1 < k < N — 2 

i—k 

n/v-i — cos &N —2 


Using the facts that 


J8{kx) 


dx 


r • 
/ si 
Jo 


sin* 9 dO — tt* —A 


k 

I W) 


and 

We have 


r(| + i) 
r(i) = nAF 


/ /* 2 tt pit rir 

8(x-n)dn = J ^ •■•J <5(|x| cos 0/v-2) sin 1 9y • • sin^ -2 9n-2^0‘ • -d9^^ 

^LiA) 

AW Ui + i)J 

(D.1) 


gN-l 


2 lx I 


_ JV-1 

2ir~r~ 


l*|r(V) 
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